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The supercooled Kob-Andersen binary Lennard-Jones (KABLJ) mixture crystaUizes partially in 
lengthy computer simulations by forming pure FCC crystals of the majority component. The general 
thermodynamic and kinetic stability of supercooled binary mixtures is discussed, emphasizing the 
importance of negative mixing enthalpy whenever present. Based on this a modification of the 
KABLJ mixture is proposed that is less prone to crystallization and faster to simulate; this is 
obtained by removing the like-particle attractions by switching to Weeks-Chandler-Andersen type 
potentials while maintaining the unlike-particle attraction. 

I. INTRODUCTION 

As computers get faster, simulations of the highly viscous liquid phase preceding glass formation become increasingly 
realistic. In this context it is nice to have a standard model system, something like the Ising model for critical 
phenomena. For several years binary binary Lennard-Jones (BLJ) mixtures have served this purpose - in particular 
the Wahnstrom and Kob-Andersen systems [1,2], because they are easy to simulate and were never found to crystallize. 
The Kob-Andersen BLJ consists of two types of Lennard-Jones particles, 80% large (A) particles and 20% small (B) 
particles. The BLJ potentials are modifications of the potentials devised by Weber and Stillinger [3], who constructed 
the pair potentials for the binary mixture based of physical-chemical data for the Ni8oP2o alloy. The Kob-Andersen 
potentials describe a strongly non-ideal mixture due to an AB attraction that is three times stronger than the BB 
attraction. This ensures a large negative mixing enthalpy (and energy) , which suppresses crystallization into pure A 
crystals as detailed below. 

Recently, both the Wahnstrom (Wa) and the Kob-Andersen (KA) systems were shown to crystallize in lengthy 
computer runs [4]. This motivated the present paper that has three purposes. First, we briefly review the general 
theory of thermodynamic and kinetic stability of supercooled binary mixtures (Sec. II). Secondly, we detail the 
crystallization of the Kob-Andersen systems (Sec. Ill), while the crystallization of the Wahnstrom liquid will be 
described elsewhere [5]. Finally, we suggest a modification of the KA binary system that is more stable and faster 
to simulate; in fact, we have not been able to crystallize the new system even in month long simulations. The idea 
is to keep the KA system's large negative mixing enthalpy, but remove the AA and BB attractions by adopting 
Weeks-Chandler- Andersen (WCA) type potentials between like particles (Sec. IV). Section V gives a brief discussion. 

II. GENERAL TREATMENT 

Consider a binary mixture of Na solvent (A) particles and Nb solute (B) particles, with total number of particles 
N = A''^ -|- Nb- The thermodynamic stability of this system against crystallization is expressed by the melting 
temperature of a pure A crystal, T£^g ^, as a function of the concentration of A particles. The latter quantity is 

conveniently expressed in terms of the fraction = Nb/N of B particles (a;^ + a^B ~ -*■)■ consider the usual 
case of externally controlled temperature T and pressure p. At the melting temperature of the pure A crystalline 
phase in the AB mixture, the chemical potential of the A particles in the crystal equals the chemical potential of 
the A particles of the liquid mixture. The chemical potential is the Gibbs free energy per particle. The change in 
Gibbs free energy per A particle at melting, AG^^g^j^g ^ can be divided into two therms, AGf^g ^ and AGj^^-^ ^, 
where AG£yg is the change in Gibbs free energy per A particle upon melting an A crystal into pure A liquid, and 
AGjjjjjj. is the change in Gibbs free energy per A particle going from pure to mixed liquid. The melting temperature 
-^fus A determined by 

^'^trans,A = ^<^fus,A(^fus,A) + ^'^mix,A(^fus,A) =0- (1) 
The case of pure A (a;g — 0) will be denoted by an asterisk, thus Tr » is the melting point of the pure A crystal 
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into pure A liquid. In this case AG^j.^^^^ ^ = AG^^^ ^ and since G = H — TS, one has AH^^^^ ^ — Tp|^^ ^ A5|^^^ ^ = 

where here and henceforth all thermodynamic quantities are per A atom. In the general case where .xg ^ 
one has AJi/f^g ^ — T£^g ^ASf^g ^ = — AGj^j^ j^. When deriving the standard expression for the freezing point 
depression one makes the approximation that crystal and liquid have the same specific heats, i.e that Ai?£^g 
and AS'£^g ^ are temperature independent. Selfconsistently with this we assume that AHj^^-^ and AS'jj^jjj. ^ 
temperature independent. Writing ^ — ^ — ATj^g ^, where ATf^g ^ is the melting point depression due 
to the presence of B particles and via AiJf^g ^ = A^'^fus A' implies that AT£^g A^'^fus A ~ ~^^mix A' 
and 



^^fus,A _ ^'^mix,A 
^fus,A ^-^fus,A 



(2) 



In the standard textbook treatment one assumes that ^H^^-^ = and that the mixing entropy is ideal, i.e., 
that the mixing entropy divided by A'' = A^i + Nb is given by AS'jj^jjj. = —kB[xj^ln{xj^) + a:gln(a;3)]. This 
expression separates into a contribution from the A particles and one from the B particles. Per A particle we thus 
have ASjjjjjj. ^ = —ks = ln(l — xg). In the limit xg ^ Eq. (2) under these assumptions reduces to 

the well-known expression 



^^fus,A ^ ^B^fus,A 
^fus.A ^^ius,A 



B- (3) 



More generally, AHjj^j^ A ~ does not apply - in fact for the Kob- Andersen liquid the mixing enthalpy is large 
(and negative) and this term cannot be ignored. The mixing entropy, on the other hand, is still realistically given by 
^•^ideal mix A- more general case Eq. (2) becomes 

^^fus,A ^ ^^mix,A + feB^fus,A H^a) 
^fus,A ^^fus,A 

A negative mixing enthalpy clearly implies a fruther melting point depression beyond that of the traditional treatments. 
If there is a large negative mixing enthalpy, this effect is important. 

The stability of a supercooled liquid mixture depends not only on its absolute (thermodynamic) stability against 
crystallization, but is also affected by kinetic effects. Thus the less supercooled the liquid is, the larger is the critical 
nucleus and the more stable is the liquid. This means that it takes longer time before a critical nucleus is generated 
by a thermal fluctuation, i.e., the supercooled liquid is more stable. At a given temperature the supercooled liquid 
is more stable the more negative AGjj^jjj. ^ becomes. To see this divide the creation of a crystal nucleus of pure A 
particles in the mixture 

7VA(mix) (crystal, a; A = 1) (5) 

into two steps: 

A'A(mix) — > A'A(liquid, a;^ = 1) 
TV^ (liquid, a; A = 1) ^ A/^ (crystal, a;^. = 1) • (6) 

The first reaction is a density fluctuation from the mixture to a pure liquid domain of A particles, the second is 

crystallization. The crystallization is (qualitatively) described by classical nucleation theory (CNT). In its simplest 
formulation the number of particles, N*, in the critical nucleus is given [6] by 
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where is the solid-hquid surface tension, pg is the crystal number density and A/x^ is the change in Gibbs free 
energy per A particle by going from the crystal to pure A liquid. When the creation of a pure A crystal takes place 

in a mixture, the lowering of Gibbs free energy by the crystallization process is reduced by 5ji = — AG^j^j^j. = 
— Ai/jjjjjj. ^ + T£^g ^ASj^jjj. ^ pa ~^A mix'^^pot ~ ^s2f^g ^ In(x^). Since the size of the critical nucleus A''* varies 
as N* oc (A/u^)~^, when A//^ is replaced by Ayit^ — 5iip^ to lowest order the change in critical nucleus size SN* is 
given by 5N*/N* = SSup^/ A^p^, thus 



SN* ^ ^A,mix«p°t(2;A) + ^srfus,A H^a) 

Ar^ = "^ Ama " ^' 

Both the negative mixing energy and the mixing entropy terms tend to increase the size of the critical nucleus. Note, 
however, that if the crystallization is instead to, e.g., an AB-crystal from an A-rich mixture with positive AB binding 
energy, the size of the critical nucleus is reduced. This prediction agrees with our simulations of crystallization in 
supercooled mixtures of varying composition. 

The above stability considerations indicate that if a binary mixtme crystallizes into, e.g., an AB-type crystal, a 
negative mixing energy will enhance the tendency of crystallization by decreasing the size of the critical nucleus. 
Confirming this, we find that increasing the fraction of B particles to .xg=0.5 for KA-type BLJ mixtures results 
in systems that more quickly crystallize into an AB (CsCl structure) crystal. This is consistent with the results of 
Fernandez and Harrowell [7] , who found that the T = equilibrium phase of the K ABL J mixture consists of coexisting 
pure A (FCC) and AB (CsCl structure) crystals. 

The mixing energy term plays also an important role for the KA mixture and here we observe (see later) that the 
nucleation time depends on details in the truncation of the pair potentials in agreement with Eq. 8. 



III. THE WAHNSTROM AND THE KOB-ANDERSEN BINARY LENNARD-JONES MIXTURES 



The simplest case of a real mixture is that given by the semi-empirical Lorentz-Berthelot (LB) rules [8] relating the 
energy parameters a and e in the pair potentials between the different species of the AB mixture: 



'^AB = (^AA + '^Bb)/2 (9) 
^AB = V^AA^BB (10) 

These rules apply for atoms with weak dispersion attractions [9], and they work well for simple mixtures of for 
instance noble-gas atoms [8], i.e., for interactions between spherically symmetric atoms with unperturbed valence- 
electron orbitals. 

There are two standard models for highly viscous supercooled liquids, the model introduced by Wahnstrom (WA) 
[1] and the model by Kob and Andersen (KA) [2]. Both models imply Lennard- Jones potentials. The WA model 
obeys the LB-rules by having 



'^AB = ('^AA + '^Bb)/2 (11) 

with CTgg = 1/1. 2(7^^ and e^g = e^A ~ ^BB ~ 1' model strongly disobeys the LB-rules: The LJ-potential 

parameters for the KA mixture are cr^g = 0-8o"aA' ''^BB ~ 0.88(7 e^g = 1-5*^AA' ^BB ~ O-Se^^. This 
results in a tendency for the small solute B particles to "glue" to the A particles. This mixture has a significant 
non-ideal mixing energy and it is considerably more stable against crystallization than the WA mixture. 

We have calculated [10] the melting point temperature in the interval £ [0, 0.2] for these two models, the result 
is shown in Figure 1. The melting point depression by mixing is traditionally given at constant pressure. Most MD 
simulations are, however performed for isochores where Tf^g ^(^^b) depends on the partial volumes of the species in 
the mixture. In the case of a KA mixture the pressure decreases with an increasing fraction of small B particles at 
constant volume and Tf-^g ^ is further depressed. 

As can be seen from the figure the negative mixing energy lowers the melting point temperature significantly. The 
Wa model is indeed more prone to crystallization than the KA model and recently Pedersen and co-workers reported 
crystallization of the Wa model after lengthy computer runs [5] . The crystal is the MgZn2 phase consisting of particles 
in the ratio 1:2, different from the 1:1 ratio defining the Wahnstrom BLJ. Thus concentration fluctuations precede (or 
at least correlate with) crystal formation. These findings inspired us to investigate whether the more widely studied 
KA model also crystallizes in lengthy computer runs. 
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FIG. 1: Melting point temperature Tjf^g of the solvent of A-particIes at the external pressure p=10 as a function of the 
particle fraction xj^ of solute particles, and in the consentration interval . xg G [0,0.2] where the system crystallizes into a 
crystal of pure (fee) A-particles. With full line is Tjf^g a(^b) ideal mixture, with dots is the Wahnstrom-model [1], with 

small dashes is for the Kob- Andersen mixture [2], and with long dashes is the present modified mixture (see Sec. IV). 




FIG. 2: (a) Particle positions projected onto a plane for the standard Kob-Andersen binary Lennard-Jones liquid (the large 
(A) particles are green, the small (B) particles are black). The system is shown after 1.5 fj,s of simulation (Argon units) at 
T=0.40 and density 1.2 (dimensionless units). After this simulation time the A particles phase separated and formed a large 
crystal, as is clear in the second figure showing the same configuration with all A particles removed. 



We performed molecular dynamics simulations of 1000 (and 10000) particles of the standard Kob-Andersen binary 
Lennard-Jones liquid in the NVE and NVT ensembles. The system was simulated at the density 1.2 in dimensionless 
units at varying temperatures. The standard time-reversible leap-frog (NVT) algorithm [11] was used with a time step 
of 0.005 (in Argon units: w 10""'^* s). The software used has been described elsewhere [12]; it utilizes a double sorting 
of neighbor particles that makes it possible to simulate 1 /is within 2-3 days of computing on a standard computer. 

For temperatures above 0.45 no crystallization was detected. In the temperature interval [0.39, 0.45] a suspicious 
drop in pressure taking place typically after w 10* time steps indicates that something is happening. It turns out that 




FIG. 3: Radial distribution functions 5aa(^) A-particles at the density p=1.2 and the temperature T=0.50 for a 

mixtures with particle fraction a;g=0.2 of B-particles. Full line is for a KA-mixture and the dashed curve is for the modified 
mixture. 




FIG. 4: The corresponding radial distribution functions b('') ^'-'^ distribution. 



the system phase separates such that the A particles form fairly large regions with no B particles present. Linked to this 
phase separation is a crystallization of the A particles. Ten simulations were performed in the [0.39, 0.45] temperature 
interval; after 4. x 10^ time steps (w 4/is) eight of these ten simulations phase separated with crystallization of the A 
particles. The general treatment of stability by supercooling a mixture indicates, 

Figure 2 shows a representative example of our results, giving the positions of the particles after crystallization at 
T=0.40 (NVT-MD). There is a large region of pure A particles showing clear crystalline order. 

that if a binary mixture can crystallize into e.g. a AB-type crystal, then the mixing energy can enhance the 
tendency of crystallization, in CNT by decreasing the size of the critical nucleus. We find that increasing the fraction 
of B-particles to a;g=0.5 results in that the KA- mixture quickly crystallizes into an AB (CsCl structure) crystal in 
accordance with the CNT prediction. The result is as mentioned also consistent with the results of Fernandez and 
Harrowell [7], who found that the T = equilibrium phase of the KA consists of a coexisting A (fee) and AB (CsCl 
structure) crystals. 

The stability osf the supercooled mixture depends on detailes in the cut-off of the potentials. Originally Kob and 
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FIG. 5: The corresponding radial distribution functions (?g b('') ^'-'^ distribution. 



Andersen used a cut-off at 2.5cra_/3 (a, (3 = A, B). If one instead of these values uses the same cut-off for all three pair 
potentials, i.e. 2.5aj^ ^, this corresponds to a shift from 2.5 aj^ g to 3.125 cr^ g. By this extension of the range of 
the AB attraction the mixing energy becomes more negative than in the original KA mixture by an amount of 0.04 
^A A' This may seem insignificant, but according to Eq. (8) the corresponding change in the critical nucleus size is of 
order 10%. Consistent with this we find that the KABLJ mixture with this cutoff' must be simulated typically three 
times longer before it crystallizes than when using the original KA parameters. 



IV. A MODIFIED BINARY MODEL SYSTEM NOT PRONE TO CRYSTALLIZATION 

With the general theory for melting point depression and nucleation in mind it is possible to modify the existing 
models in such a way that they are not prone to crystallization. We want to maintain the success of the existing 
models, which are set up to simulate real supercooled mixtures, in the case of KA-model a binary mixture of Ni8oP2o 
alloy. For this reason the sizes of the molecules in the KA-mixture, given by the an values, are maintained. We also 
want to have the same structure of the KA-mixture, given by their correlation functions. A modified KA-mixture 
with only repulsive forces between A-A and B-B particles, but where the attraction between A- and B-particles is 
maintained 



with 



Uij {Vij ) 



4e,: 



(^) -(^) 



ui^jinjicut)), r^j<rij(cut) 
r > rij{cut), 



(12) 



ri,i(cut) = 2i/Vi,i 
ta^bM) = 2.5aA,B 



(13) 
(14) 



meets these requirements. A priory it is to be expected [13] that the structure of the mixture only deviates marginally 

from the originally KA-model because the repulsive parts of the potential interactions are not changed from the original 
KA-model. The modified binary model is indeed much more stable and we have so far not been able to crystallizate 
the mixture (for a particle fraction .Tg = 0.2). There are two reasons for this resistance against crystallization by 
supercooling. Firstly the melting point temperature, ^, of the pure solvent particles is lowered significantly (see 

Figure 1) and secondly the relative importance of the mixing energy, mix^pot(^B) ' bigger due to a stronger 
violation of the LB-rules for the energy and thereby a lower melting point temperature and an increase of the size 
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FIG. 6: log — log plot of the mean square displacement for A-particles as a function of time (in unit a^^m/e « 2. x 10 '^^s)and 
for different temperatures (from the left): T=1.Q0, 0.40, 0.35, 0.325, 0.30, 0.275 and 0.25. respectively. 



of the critical nucleus. Another advantage of the modified system is that it is much faster to simulate by using a 
two-step sorting of nearest neighbor [12]. In additional to that one can use a bigger value of the time increment in 
the MD simulation due to the lower temperature, so not only does the mixture not crystallizes; but it can also be 
cooled futher down under equilibrium condition (see Figure 6-8) and followed a much longer times (0.1ms). 

The structure of the modified mixture is, as expected, almost the same as in the KA-mixture. Figure 3-5 show 
the radial distribution functions, ga,/3(r) for the modified mixture and for the KA-mixture. The similarity is well 
known. (The so-called WCA- system [13] with only repulsive LJ forces has been used in perturbation theories for 
dense liquid.) As can be seen from the figures only the distribution, b('^) '^-'^ ^^^^ B-particles is affected. This 
change in g-Q ^^{r) is due to the relative increased "binding-energy" between the solvent A-particles and the solute 
of B-particles which, however, is not so visible in g (r) due to the many A-B configurations (notice the diflierent 
y-scales of the two figures). 

The KA- and the modified mixture were cooled down from T=1.5. The KA- mixture crystallizes as mentioned 
in the temperature interval T 6 [0.39,0.45]; but it is possible to quench the KA-system down to T=0.375 and still 
determine the self diffusion constant, D. The modified mixture can be cooled down to T=0.25, where the system was 
followed 10^° « one 0.1ms time step and with a diffusion constant for the A-particles, D(A) w 1. x 10^°**. The B- 
particles behaves in a similar way; but as the temperature is decreased the ratio D(B)/D(A) changes from 1.7 for 
T=1.5 to 9.7 for T=0.25. The values of the self-diffusion constants D(T) were obtained from the slopes of the mean 
square displacements as a function of the time. A log log-log plot of D(T) for different temperatures are shown in 
Figure 6. At the low temperatures the ballistic- and diffusive time-regimes are very well separated. 

The diffusion const ant.s for A- and B-particles are shown in Figure 7 together with low temperature data for the 
KA-mixture (A-particles) [14]. It is possible to scale the diffusion constants for the two mixtures as is demonstrated in 
Figure 8, where log(D) for the KA-mixture (-|-) connected with dots, is plotted as a function of 1/(T*0.81). A similar 
behavior for "fragility invariance" for systems with the same repulsive potentials has been observed before [15]. 

V. DISCUSSION 

The highly viscous fluid state below the equilibrium melting point temperature is found in many systems ranging 
from simple organic one-component systems and alloys to complex biochemical substances and ionic liquids. Along 
with the increasingly interest in this state of matter it is important to have models not prone to crystallisation. The 
general thermodynamics theory for (equilibrium) melting and classical nucleation theory gives an indication of how 
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FIG. 7: An Arrhenius plot. logD{l/T). of the self diffusion constant D. With + is D(A) for the KA mixture, and the points 
given by X and connected with lines is D(A) for the modified binary mixture (see III B). The diffusion constants, D(B) for the 
smaller B-particles in the modified binary mixture are shown with 
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FIG. 8: With full line is the logT){l/T) for the modified mixture ( the data points from Fig.8 are just connected by straight 
lines) and with dashes are the corresponding scaled KA-data log(D){^j-^). 



to create models which suppress crystallisation by cooling. We have given one simple example by modifying the 
Kob- Andersen model which is the most used MD model for binary highly viscous fluidis but it is straight forward to 
extend the model to include other interactions and a more refined energy functions without loose its stability against 
crystallization. A general recipe of energy functions for stable highly viscous liquids are, however difficult to set up 
since many mixtures exhibits "reacting systems" with crystals with mixed compositions where an increased binding-or 
mixing energy can enhance crystallization of e.g. A,B crystals. Both the Kob-Andersen and the Wahnstrom models 
are examples of this phenomenon. 

The modified model enhanches the relative attraction between unlike species by removing the attraction between 
the A,A- and B,B-interactions. Supercooled liquid (WCA-) mixtures with solely repulsive LJ-potentials have been 
reported before [16], [17]; but we observe that both systems [16], [17] phase separate and crystallized by cooling. 
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By maintaining the attraction between unlike species the highly viscous state is stabilized. The modified mixture 
disobeys the LB-energy rule much more than the KA-mixture. The LB-energy rule is a geometric mean of interaction 

energy between two particles with energy scaling parameters, ep^^ and egg. The two models have the same energy 
scaling parameters, but a much stronger violation of the energy rule in the interaction interval 6 oo] 
where the WCA-potentials are zero. The result of this "cut and shift" is that the B-particles are almost "covalent" 
bounded to the A-particles and the violation of the LB-rule for the mixing energy gives recipe of creating fragile and 
stable supercooled liquid mixtures. 
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